global pollutants "ethanol_ozone_data"
global fuelprices "fuel_prices_data"

set matsize 10000


************************************************************************************
************************************************************************************
************************************************************************************
* This program should be in memory since it is called on when opening ethanol_ozone_data.dta
capture program drop prepareVariables
program define prepareVariables
************************************************************************************
************************************************************************************
************************************************************************************

* Generate logarithms of pollutant concentrations
gen logO3=log(O3)
* When concentration for other pollutants is 0, instead of missing take as the logarithm the lowest value within site-discontinuity
foreach var1 in PM25 CO NO NO2 NOX {
	gen log`var1'=log(`var1')
	bysort dis siteid: egen log`var1'_min=min(log`var1') if hour>=7&hour<=11
	replace log`var1'=log`var1'_min if `var1'==0&hour>=7&hour<=11 
	drop log`var1'_min
	bysort dis siteid: egen log`var1'_min=min(log`var1') if hour>=17&hour<=20
	replace log`var1'=log`var1'_min if `var1'==0&hour>=17&hour<=20
	drop log`var1'_min
	}

* Generate logarithms of meteorological conditions
foreach var1 in tpl0 rdl0 hml0 wsl0 {
	gen log`var1'=log(`var1')
	}

* Generate mean precipitation in the contemporaneous hour to three hours earlier, in mm/h.
bysort dis siteid (date hour): gen ppl0to3=(ppl0+ppl0[_n-1]+ppl0[_n-2]+ppl0[_n-3])/4

* Generate site-date level mean road traffic congestion during the morning commuting hours and during the evening commuting hours
gen morning_traffic=kmcityl0 if hour>=7&hour<=11
bysort dis siteid date: egen kmcity_am=mean(morning_traffic)
gen logkmcity_am=log(kmcity_am+1)
gen afternoon_traffic=kmcityl0 if hour>=17&hour<=20
bysort dis siteid date: egen kmcity_pm=mean(afternoon_traffic)
gen logkmcity_pm=log(kmcity_pm+1)
drop morning_traffic afternoon_traffic

* Generate maximum 8-hour average for ozone (always missing for hour 6 due to instrument calibration)
sort dis siteid time hour
replace O3=0 if hour==6
gen mean_O3=(O3+O3[_n+1]+O3[_n+2]+O3[_n+3]+O3[_n+4] ///
		+O3[_n+5]+O3[_n+6]+O3[_n+7]+O3[_n+8])/8 if hour<6
replace mean_O3=(O3+O3[_n+1]+O3[_n+2]+O3[_n+3]+O3[_n+4] ///
		+O3[_n+5]+O3[_n+6]+O3[_n+7])/8 if hour<=16 & hour>=7
replace O3=. if hour==6
bysort dis siteid date: egen max_mean_O3=max(mean_O3)
gen hour_max_start=hour if (max_mean_O3-mean_O3<=0.0000001)
bysort dis siteid date: egen period_start=max(hour_max_start)
gen period_end=period_start+8 if period_start<=6
replace period_end=period_start+7 if period_start>6
drop mean_O3 hour_max_start

end


************************************************************************************
************************************************************************************
************************************************************************************
capture program drop summaryStatistics
program define summaryStatistics
************************************************************************************
************************************************************************************
************************************************************************************

use $pollutants, clear
* Pollutant concentrations
sum O3 if hour>=7&hour<=11
sum O3 if hour>=12&hour<=16
sum e25 if hour>=12&hour<=16&O3~=.
sum PM25 if hour>=7&hour<=11
sum PM25 if hour>=17&hour<=20
sum CO if hour>=7&hour<=11
sum CO if hour>=17&hour<=20
sum NOX if hour>=7&hour<=11		// Table 2 includes siteid 34 which was excluded from regression
sum NOX if hour>=17&hour<=20	// due to not meeting 70% availability criteria
* Meteorology
preserve
	collapse (mean) tpl0 rdl0 hml0 wsl0 ppl0, by (date hour dis)
	foreach var1 in tpl0 rdl0 hml0 wsl0 ppl0 {
		sum `var1' if hour>=7 & hour<=11
		sum `var1' if hour>=12 & hour<=16
	}
restore
* Wind direction: Four weather stations measured wind direction already in 2010
preserve
	keep if ((siteid==2|siteid==5|siteid==17|siteid==27)&(dis~=4)) | ///
			((siteid==2|siteid==5|siteid==17           )&(dis==4))
	sum wdqdt1l0 wdqdt2l0 wdqdt3l0 wdqdt4l0 if hour>=7&hour<=11
	sum wdqdt1l0 wdqdt2l0 wdqdt3l0 wdqdt4l0 if hour>=12&hour<=16		
restore
* Wind direction: More weather stations measured wind direction by 2013
preserve 
	keep if ((siteid==2|siteid==5|                     siteid==17|           siteid==27                                 )&(dis~=4)) | ///
			((siteid==2|siteid==5|siteid==7|siteid==16|siteid==17|siteid==28|           siteid==35|siteid==36|siteid==37)&(dis==4))	
	sum wdqdt1l0 wdqdt2l0 wdqdt3l0 wdqdt4l0 if hour>=7&hour<=11
	sum wdqdt1l0 wdqdt2l0 wdqdt3l0 wdqdt4l0 if hour>=12&hour<=16		
restore
* Thermal inversion
preserve
	collapse (mean) dvti0to199 dvti200to499 if hour==9|hour==21, by(date hour dis)
	sum dvti0to199 dvti200to499 if hour==9
	sum dvti0to199 dvti200to499 if hour==21
restore
* Road congestion
preserve
	collapse (mean) kmcityl0, by(date hour dis)
	sum kmcityl0 if hour>=7&hour<=11
	sum kmcityl0 if hour>=17&hour<=20
restore
* Ethanol-to-gasoline price ratio
preserve
	collapse (mean) pe_pg_ratio, by(date dis)
	sum  pe_pg_ratio
restore

end


************************************************************************************
************************************************************************************
************************************************************************************
capture program drop tables
program define tables
************************************************************************************
************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 3
use $pollutants, clear
prepareVariables
local sitesforregdis1 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis2 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis3 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis4 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15 | siteid==16 | siteid==18 | siteid==22 | siteid==27 | siteid==28 | siteid==29 | siteid==31 | siteid==33 | siteid==35 | siteid==37"

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==1&(`sitesforregdis1'), cluster(sitedate)
dis "First discontinuity" ", " "Coef:" -_b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==2&(`sitesforregdis2'), cluster(sitedate)
dis "Second discontinuity" ", " "Coef:" _b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==3&(`sitesforregdis3'), cluster(sitedate)
dis "Third discontinuity" ", " "Coef:" -_b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==4&(`sitesforregdis4'), cluster(sitedate)
dis "Fourth discontinuity" ", " "Coef:" _b[D] ", " "SE:" _se[D]

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 4A
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
preserve
quietly{
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(30) rho(0.677) all /*CCT*/
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(31) rho(0.853) all /*IK*/
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(25) rho(0.691) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(35) rho(0.910) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(22) rho(0.610) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(24) rho(0.598) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(28) rho(0.512) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(23) rho(0.582) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
}

* Results for Table 4A
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st  ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd  ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd  ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th  ", SE: " se_ik_4th
scalar drop _all
restore

/* Below are original regressions with CCT and IK criteria. We consider integer days.
We chop off the decimal with decimal value smaller than 0.7 and add one more day otherwise.
For the remaining sections, we use the same method to select bandwidth. */
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, bwselect(CCT) all 
rdrobust res_o3_1st time if dis==1, bwselect(IK) all 

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 4B
use $pollutants, clear
prepareVariables
quietly{
gen res_o3max_1st=.
gen res_o3max_2nd=.
gen res_o3max_3rd=.
gen res_o3max_4th=.
bysort dis siteid date: egen max_O3=max(O3) if hour>=12&hour<=16
keep if abs(O3-max_O3)<0.00001
gen logO3_max=log(max_O3)

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
quietly reg logO3_max c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3max_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
quietly reg logO3_max c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3max_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
quietly reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3max_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
quietly reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3max_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
rdrobust res_o3max_1st time if dis==1, h(27) all
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3max_1st time if dis==1, h(32) all
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3max_2nd time if dis==2, h(27) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional] 
rdrobust res_o3max_2nd time if dis==2, h(31) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional] 

rdrobust res_o3max_3rd time if dis==3, h(26) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3max_3rd time if dis==3, h(26) all 
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3max_4th time if dis==4, h(29) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3max_4th time if dis==4, h(28) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
}

* Results for Table 4B
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st  ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd  ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd  ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th  ", SE: " se_ik_4th
scalar drop _all

* CCT and IK bandwidth criteria
rdrobust res_o3max_1st time if dis==1, bwselect(CCT) all 
rdrobust res_o3max_1st time if dis==1, bwselect(IK) all

rdrobust res_o3max_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3max_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3max_3rd time if dis==3, bwselect(CCT) all 
rdrobust res_o3max_3rd time if dis==3, bwselect(IK) all

rdrobust res_o3max_4th time if dis==4, bwselect(CCT) all 
rdrobust res_o3max_4th time if dis==4, bwselect(IK) all

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 5
use $pollutants, clear
prepareVariables
quietly{
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
gen res_o3=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio )#i.hour ///
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
preserve 
quietly{
collapse res time_1, by(date siteid dis)
rdrobust res time_1, h(23) all /*CCT*/
scalar b_cct=_b[Conventional]
scalar se_cct=_se[Conventional]
rdrobust res time_1, h(27) all /*IK*/
scalar b_ik=_b[Conventional]
scalar se_ik=_se[Conventional]
}

* Results for Table 5
dis "CCT Criterion: " "Coef: " b_cct ", SE: " se_cct
dis "IK Criterion: "  "Coef: " b_ik  ", SE: " se_ik
restore

* CCT and IK bandwidth criteria
preserve 
collapse res time_1, by(date siteid dis)
rdrobust res time_1, bwselect(CCT) all
rdrobust res time_1, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 6A
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=7&hour<=11&dis==1
predict res_temp if (siteid==`var1')&hour>=7&hour<=11&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=7&hour<=11&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=7&hour<=11&dis==2
predict res_temp if (siteid==`var1')&hour>=7&hour<=11&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=7&hour<=11&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=7&hour<=11&dis==3				  
predict res_temp if (siteid==`var1')&hour>=7&hour<=11&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=7&hour<=11&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=7&hour<=11&dis==4	
predict res_temp if (siteid==`var1')&hour>=7&hour<=11&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=7&hour<=11&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(27) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(33) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(24) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(25) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(34) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(32) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(36) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(29) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(25) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(28) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results for Table 6A
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_o3* time, by(date siteid dis)

rdrobust res_o3_1st time if dis==1, bwselect(CCT) all 
rdrobust res_o3_1st time if dis==1, bwselect(IK) all 

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 6B
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=11&hour<=18&dis==1
predict res_temp if (siteid==`var1')&hour>=11&hour<=18&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=11&hour<=18&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=11&hour<=18&dis==2
predict res_temp if (siteid==`var1')&hour>=11&hour<=18&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=11&hour<=18&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=11&hour<=18&dis==3				  
predict res_temp if (siteid==`var1')&hour>=11&hour<=18&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=11&hour<=18&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=11&hour<=18&dis==4	
predict res_temp if (siteid==`var1')&hour>=11&hour<=18&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=11&hour<=18&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(30) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(33) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(24) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(32) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(21) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(32) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(28) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(27) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(24) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(27) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results for Table 6B
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_o3* time, by(date siteid dis)

rdrobust res_o3_1st time if dis==1, bwselect(CCT) all 
rdrobust res_o3_1st time if dis==1, bwselect(IK) all 

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 6C
/*
Results differ slightly from what we present in Table 6D due to a minor
revision in how we generate the maximum 8-average for ozone (see code above).
*/
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
			     c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==1
predict res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==2
predict res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==3				  
predict res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==4	
predict res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(30) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(33) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(22) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(34) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(26) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(23) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(31) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(23) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(25) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(25) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}


* Results for Table 6C
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_o3* time, by(date siteid dis)

rdrobust res_o3_1st time if dis==1, bwselect(CCT) all 
rdrobust res_o3_1st time if dis==1, bwselect(IK) all 

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 6D
/*
Results differ slightly from what we present in Table 6D due to a minor
revision in how we generate the maximum 8-average for ozone (see code above).
*/
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg O3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
			     c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==1
predict res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg O3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==2
predict res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg O3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==3				  
predict res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg O3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==4	
predict res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=period_start & hour<=period_end&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(23) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(24) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(36) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(42) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(26) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(26) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(26) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(30) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(27) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(28) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results for Table 6D
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_o3* time, by(date siteid dis)

rdrobust res_o3_1st time if dis==1, bwselect(CCT) all 
rdrobust res_o3_1st time if dis==1, bwselect(IK) all 

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 7A
use $pollutants, clear
prepareVariables
quietly{
gen res_pm25_3rd=.
gen res_pm25_4th=.

foreach var1 in 8 {
quietly reg logPM25 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==3
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==3, residual
replace res_pm25_3rd=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==3
drop res_temp
}
foreach var1 in 8 27 31 {
quietly reg logPM25 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==4
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==4, residual
replace res_pm25_4th=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_pm25* time, by(date siteid dis)
rdrobust res_pm25_3rd time if dis==3, h(25) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_pm25_3rd time if dis==3, h(52) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]
rdrobust res_pm25_4th time if dis==4, h(29) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_pm25_4th time if dis==4, h(48) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_pm25_3rd if dis==3
replace res_all=res_pm25_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(33) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(48) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results for Table 7A
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_pm25* time, by(date siteid dis)
rdrobust res_pm25_3rd time if dis==3, bwselect(CCT) all
rdrobust res_pm25_3rd time if dis==3, bwselect(IK) all

rdrobust res_pm25_4th time if dis==4, bwselect(CCT) all
rdrobust res_pm25_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_pm25_3rd if dis==3
replace res_all=res_pm25_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 7B
/*
The results below for the second and third discontinuities differ slightly from
those presented in Table 7B. When preparing Table 7B, site 32 was not included
in the discontinuity #2 sample, and site 17 was not included in the discontinuity
#3 sample, due to data not being available at some point in the project. The data have
now been supplemented. To replicate Table 7B, including the estimates for the pooled
sample, drop these sites from the respective estimation samples. The bandwidth should  
be changed to that displayed in Table 7B.
*/
use $pollutants, clear
prepareVariables
quietly{
gen res_co_1st=.
gen res_co_2nd=.
gen res_co_3rd=.
gen res_co_4th=.

foreach var1 in 1 3 5 7 8 10 17 20 27 29 31 32 {
reg logCO (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==1
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==1, residual
replace res_co_1st=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==1
drop res_temp
}
foreach var1 in 1 3 5 7 8 10 17 20 27 29 31 32{
reg logCO (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==2
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==2, residual
replace res_co_2nd=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==2
drop res_temp
}
foreach var1 in 1 3 5 7 8 10 17 20 27 29 31 {
reg logCO (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==3
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==3, residual
replace res_co_3rd=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==3
drop res_temp
}
foreach var1 in 1 3 5 7 8 10 16 17 20 27 28 29 31 32 36 {
reg logCO (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==4
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==4, residual
replace res_co_4th=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_co* time, by(date siteid dis)
rdrobust res_co_1st time if dis==1, h(22) all
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_co_1st time if dis==1, h(90) all
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_co_2nd time if dis==2, h(39) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_co_2nd time if dis==2, h(38) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_co_3rd time if dis==3, h(29) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_co_3rd time if dis==3, h(31) all 
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_co_4th time if dis==4, h(26) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_co_4th time if dis==4, h(27) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_co_1st if dis==1
replace res_all=res_co_2nd if dis==2
replace res_all=res_co_3rd if dis==3
replace res_all=res_co_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(33) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(28) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results for Table 7B
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_co* time, by(date siteid dis)
rdrobust res_co_1st time if dis==1, bwselect(CCT) all
rdrobust res_co_1st time if dis==1, bwselect(IK) all

rdrobust res_co_2nd time if dis==2, bwselect(CCT) all
rdrobust res_co_2nd time if dis==2, bwselect(IK) all

rdrobust res_co_3rd time if dis==3, bwselect(CCT) all 
rdrobust res_co_3rd time if dis==3, bwselect(IK) all 

rdrobust res_co_4th time if dis==4, bwselect(CCT) all
rdrobust res_co_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_co_1st if dis==1
replace res_all=res_co_2nd if dis==2
replace res_all=res_co_3rd if dis==3
replace res_all=res_co_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table 7C
use $pollutants, clear
prepareVariables
quietly{
gen res_noX_1st=.
gen res_noX_2nd=.
gen res_noX_3rd=.
gen res_noX_4th=.

foreach var1 in 1 5 7 8 10 17 20 22 27 29 31 {
quietly reg logNOX (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499  ///
				  c.logkmcity_pm#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==1
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==1, residual
replace res_noX_1st=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==1
drop res_temp
}
foreach var1 in 1 5 7 8 10 17 20 22 27 29 31 {
quietly reg logNOX (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==2
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==2, residual
replace res_noX_2nd=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==2
drop res_temp
}
foreach var1 in 1 5 7 8 10 17 20 22 27 29 31 {
quietly reg logNOX (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==3
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==3, residual
replace res_noX_3rd=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==3
drop res_temp
}
foreach var1 in  1 5 7 10 17 20 22 27 28 29 31 35 36 37 {
quietly reg logNOX (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_pm#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=17&hour<=20&dis==4
predict res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==4, residual
replace res_noX_4th=res_temp if (siteid==`var1')&hour>=17&hour<=20&dis==4
drop res_temp
}
}
preserve
collapse res_noX* time, by(date siteid dis)
rdrobust res_noX_1st time if dis==1, h(24) all
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_noX_1st time if dis==1, h(51) all
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_noX_2nd time if dis==2, h(32) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_noX_2nd time if dis==2, h(53) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_noX_3rd time if dis==3, h(20) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_noX_3rd time if dis==3, h(41) all 
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_noX_4th time if dis==4, h(35) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_noX_4th time if dis==4, h(47) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_noX_1st if dis==1
replace res_all=res_noX_2nd if dis==2
replace res_all=res_noX_3rd if dis==3
replace res_all=res_noX_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(22) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(47) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore

* Results for Table 7C
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_noX* time, by(date siteid dis)
rdrobust res_noX_1st time if dis==1, bwselect(CCT) all
rdrobust res_noX_1st time if dis==1, bwselect(IK) all

rdrobust res_noX_2nd time if dis==2, bwselect(CCT) all 
rdrobust res_noX_2nd time if dis==2, bwselect(IK) all 

rdrobust res_noX_3rd time if dis==3, bwselect(CCT) all 
rdrobust res_noX_3rd time if dis==3, bwselect(IK) all 

rdrobust res_noX_4th time if dis==4, bwselect(CCT) all
rdrobust res_noX_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_noX_1st if dis==1
replace res_all=res_noX_2nd if dis==2
replace res_all=res_noX_3rd if dis==3
replace res_all=res_noX_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table A.1
* See above for 8th order polynomial
* 9th order polynomial
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.
gen time_9=time^9

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time c.time_9)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time c.time_9)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time c.time_9)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time c.time_9)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(30) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(31) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(25) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(35) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(22) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(24) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(28) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(23) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(23) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(27) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results with 9th order polynomial
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* 7th order polynomial
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(30) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(31) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(25) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(35) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(22) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(24) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(28) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(23) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(23) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(27) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results with 7th order polynomial
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* 6th order polynomial
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(30) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(31) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(25) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(35) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(22) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(24) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(28) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(23) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(23) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(27) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results with 6th order polynomial
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* No polynomial 
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(30) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(31) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(25) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(35) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(22) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(24) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(28) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(23) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(23) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(27) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}
* Results with no polynomial
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

************************************************************************************
************************************************************************************

* This part reproduces the standard errors shown in Table A.2
* This part reproduces results for the first discontinuity when the bandwidth is re-selected
* according to each bootstrap sample using the CCT criterion (thus the bootstrap sample
* consists of 12 sites * 180 dates = 2160 site-date clusters). To complete the procedure,
* compute the standard deviation of _b[Conventional] which is recorded in the log file.
log using test_bs_cct_1st, text
* The following loops 400 times, once for each bootstrap sample
forvalues var1=28769354(11)28773743{
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
keep if siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 ///
		| siteid==7 | siteid==15 | siteid==18 | siteid==22 | siteid==27 ///
		| siteid==29 | siteid==31
keep if dis==1
set seed `var1'
bsample 2160, cluster(sitedate)
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
collapse res time, by(date siteid)
rdrobust res time, bwselect(CCT) all
}
dis _b[Conventional]
}
log close

* This part reproduces results for the first discontinuity when the bandwidth is fixed
* at the Table 4A value using the CCT criterion, based on the original sample.
log using test_bs_fixedcct_1st, text
forvalues var1=28769354(11)28773743{
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
keep if siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 ///
		| siteid==7 | siteid==15 | siteid==18 | siteid==22 | siteid==27 ///
		| siteid==29 | siteid==31
keep if dis==1
set seed `var1'
bsample 2160, cluster(sitedate)
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
collapse res time, by(date siteid)
rdrobust res time, h(30) all
}
dis _b[Conventional]
}
log close

************************************************************************************
************************************************************************************

*This part reproduces the results shown in Table A.3
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
			     c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(29) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(31) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(24) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(33) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(26) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(31) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(27) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(25) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore
}

* Results for Table A.3
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st  ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd  ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd  ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th  ", SE: " se_ik_4th
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, bwselect(CCT) all
rdrobust res_o3_1st time if dis==1, bwselect(IK) all

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all 
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all 

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table A.4
* Linear phase-in over three days
use $pollutants, clear
prepareVariables
local sitesforregdis1 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis2 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis3 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis4 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15 | siteid==16 | siteid==18 | siteid==22 | siteid==27 | siteid==28 | siteid==29 | siteid==31 | siteid==33 | siteid==35 | siteid==37"
replace D=0.33 if time==0
replace D=0.67 if time==1
quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==1&(`sitesforregdis1'), cluster(sitedate)
dis "First discontinuity" ", " "Coef:" -_b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==2&(`sitesforregdis2'), cluster(sitedate)
dis "Second discontinuity" ", " "Coef:" _b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==3&(`sitesforregdis3'), cluster(sitedate)
dis "Third discontinuity" ", " "Coef:" -_b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==4&(`sitesforregdis4'), cluster(sitedate)
dis "Fourth discontinuity" ", " "Coef:" _b[D] ", " "SE:" _se[D]

* Linear phase-in over five days
use $pollutants, clear
prepareVariables
local sitesforregdis1 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis2 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis3 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15              | siteid==18 | siteid==22 | siteid==27              | siteid==29 | siteid==31"
local sitesforregdis4 "siteid==1 | siteid==2 | siteid==3 | siteid==5 | siteid==6 | siteid==7 | siteid==15 | siteid==16 | siteid==18 | siteid==22 | siteid==27 | siteid==28 | siteid==29 | siteid==31 | siteid==33 | siteid==35 | siteid==37"
replace D=0.2 if time==0
replace D=0.4 if time==1
replace D=0.6 if time==2
replace D=0.8 if time==3
quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==1&(`sitesforregdis1'), cluster(sitedate)
dis "First discontinuity" ", " "Coef:" -_b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==2&(`sitesforregdis2'), cluster(sitedate)
dis "Second discontinuity" ", " "Coef:" _b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==3&(`sitesforregdis3'), cluster(sitedate)
dis "Third discontinuity" ", " "Coef:" -_b[D] ", " "SE:" _se[D]

quietly reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.sitehour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.sitehour ///
				  i.dayofweek##i.sitehour ///
				  D if hour>=12&hour<=16&dis==4&(`sitesforregdis4'), cluster(sitedate)
dis "Fourth discontinuity" ", " "Coef:" _b[D] ", " "SE:" _se[D]

************************************************************************************
************************************************************************************

* This part reproduces the results shown in Table A.5
* Meteorology in levels
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.tpl0##c.tpl0 c.rdl0##c.rdl0 c.hml0##c.hml0 c.wsl0##c.wsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.tpl0##c.tpl0 c.rdl0##c.rdl0 c.hml0##c.hml0 c.wsl0##c.wsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.tpl0##c.tpl0 c.rdl0##c.rdl0 c.hml0##c.hml0 c.wsl0##c.wsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.tpl0##c.tpl0 c.rdl0##c.rdl0 c.hml0##c.hml0 c.wsl0##c.wsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(30) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(32) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(24) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(33) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(27) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(24) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(27) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(19) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(23) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(25) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results for meteorology in levels
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, bwselect(CCT) all
rdrobust res_o3_1st time if dis==1, bwselect(IK) all

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all 
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all 

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

* Discretize temperature: 1 C bins
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.
egen tpl0cat2 = cut(tpl0), at(20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36)

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (i.tpl0cat2 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (i.tpl0cat2 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (i.tpl0cat2 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (i.tpl0cat2 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(31) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(33) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(30) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(45) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(26) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(22) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(19) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(32) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(21) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(27) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results for temperature bins
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, bwselect(CCT) all
rdrobust res_o3_1st time if dis==1, bwselect(IK) all

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all 

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore

* Discretize radiation: 150 W/m^2 bins
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.
egen rdl0cat2 = cut(rdl0), at(150 300 450 600 750 900 1150)

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 i.rdl0cat2 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 i.rdl0cat2 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 i.rdl0cat2 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 i.rdl0cat2 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
quietly{
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, h(33) all 
scalar b_cct_1st=-_b[Conventional]
scalar se_cct_1st=_se[Conventional]
rdrobust res_o3_1st time if dis==1, h(33) all 
scalar b_ik_1st=-_b[Conventional]
scalar se_ik_1st=_se[Conventional]

rdrobust res_o3_2nd time if dis==2, h(28) all
scalar b_cct_2nd=_b[Conventional]
scalar se_cct_2nd=_se[Conventional]
rdrobust res_o3_2nd time if dis==2, h(35) all
scalar b_ik_2nd=_b[Conventional]
scalar se_ik_2nd=_se[Conventional]

rdrobust res_o3_3rd time if dis==3, h(20) all
scalar b_cct_3rd=-_b[Conventional]
scalar se_cct_3rd=_se[Conventional]
rdrobust res_o3_3rd time if dis==3, h(24) all
scalar b_ik_3rd=-_b[Conventional]
scalar se_ik_3rd=_se[Conventional]

rdrobust res_o3_4th time if dis==4, h(18) all
scalar b_cct_4th=_b[Conventional]
scalar se_cct_4th=_se[Conventional]
rdrobust res_o3_4th time if dis==4, h(21) all
scalar b_ik_4th=_b[Conventional]
scalar se_ik_4th=_se[Conventional]
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, h(22) all
scalar b_cct_pool=_b[Conventional]
scalar se_cct_pool=_se[Conventional]
rdrobust res_all time_1, h(20) all
scalar b_ik_pool=_b[Conventional]
scalar se_ik_pool=_se[Conventional]
restore
}

* Results with radiation bins
dis "First dis CCT Criterion: "  "Coef: " b_cct_1st  ", SE: " se_cct_1st
dis "First dis IK Criterion: "   "Coef: " b_ik_1st   ", SE: " se_ik_1st
dis "Second dis CCT Criterion: " "Coef: " b_cct_2nd  ", SE: " se_cct_2nd
dis "Second dis IK Criterion: "  "Coef: " b_ik_2nd   ", SE: " se_ik_2nd
dis "Third dis CCT Criterion: "  "Coef: " b_cct_3rd  ", SE: " se_cct_3rd
dis "Third dis IK Criterion: "   "Coef: " b_ik_3rd   ", SE: " se_ik_3rd
dis "Fourth dis CCT Criterion: " "Coef: " b_cct_4th  ", SE: " se_cct_4th
dis "Fourth dis IK Criterion: "  "Coef: " b_ik_4th   ", SE: " se_ik_4th
dis "Pooled CCT Criterion: "     "Coef: " b_cct_pool ", SE: " se_cct_pool
dis "Pooled IK Criterion: "      "Coef: " b_ik_pool  ", SE: " se_ik_pool
scalar drop _all

* CCT and IK bandwidth criteria
preserve
collapse res_o3* time, by(date siteid dis)
rdrobust res_o3_1st time if dis==1, bwselect(CCT) all
rdrobust res_o3_1st time if dis==1, bwselect(IK) all

rdrobust res_o3_2nd time if dis==2, bwselect(CCT) all
rdrobust res_o3_2nd time if dis==2, bwselect(IK) all

rdrobust res_o3_3rd time if dis==3, bwselect(CCT) all 
rdrobust res_o3_3rd time if dis==3, bwselect(IK) all

rdrobust res_o3_4th time if dis==4, bwselect(CCT) all
rdrobust res_o3_4th time if dis==4, bwselect(IK) all
restore

preserve
gen res_all=res_o3_1st if dis==1
replace res_all=res_o3_2nd if dis==2
replace res_all=res_o3_3rd if dis==3
replace res_all=res_o3_4th if dis==4
gen time_1=time
replace time_1=-1-time if dis==1 | dis==3
collapse res_all time_1, by(date siteid dis)
rdrobust res_all time_1, bwselect(CCT) all
rdrobust res_all time_1, bwselect(IK) all
restore


***********************************************************************************
***********************************************************************************

* This part reproduces the results shown in Table A.6 
use $fuelprices, clear
sort brand
gen brand_index=_n if brand~=brand[_n-1]
replace brand_index=brand_index[_n-1] if brand==brand[_n-1]
sort date_surveyed
gen running=date_surveyed-date("01-Feb-2010", "DMY")
keep if running>=-90&running<=89
drop if zip3==. 
quietly reg pg c.running##c.running##c.running##c.running##c.running##c.running##c.running##c.running ///
	   i.brand_index i.zip3
predict pg_resid, residual
rdrobust pg_resid running, h(40) all
scalar bpg_cct_1st_1=-_b[Conventional]
scalar sepg_cct_1st_1=_se[Conventional]
rdrobust pg_resid running, h(81) all
scalar bpg_ik_1st_1=-_b[Conventional]
scalar sepg_ik_1st_1=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pg_resid running, bwselect(CCT) all
rdrobust pg_resid running, bwselect(IK) all
preserve
gen zip2=(zip3-mod(zip3,10))/10
collapse pg_resid, by(running zip2 brand_index)
rdrobust pg_resid running, h(42) all
scalar bpg_cct_1st_2=-_b[Conventional]
scalar sepg_cct_1st_2=_se[Conventional]
rdrobust pg_resid running, h(60) all
scalar bpg_ik_1st_2=-_b[Conventional]
scalar sepg_ik_1st_2=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pg_resid running, bwselect(CCT) all
rdrobust pg_resid running, bwselect(IK) all
restore
* Results for Table A.6, price of E20/E25 gasoline, first discontinuity
dis "Coef: " bpg_cct_1st_1 ", SE: " sepg_cct_1st_1 
dis "Coef: " bpg_cct_1st_2 ", SE: " sepg_cct_1st_2 
dis "Coef: " bpg_ik_1st_1  ", SE: " sepg_ik_1st_1 
dis "Coef: " bpg_ik_1st_2  ", SE: " sepg_ik_1st_2 

use $fuelprices, clear
sort brand
gen brand_index=_n if brand~=brand[_n-1]
replace brand_index=brand_index[_n-1] if brand==brand[_n-1]
sort date_surveyed
gen running=date_surveyed-date("01-May-2010", "DMY")
keep if running>=-90&running<=89
drop if zip3==. 
quietly reg pg c.running##c.running##c.running##c.running##c.running##c.running##c.running##c.running ///
	   i.brand_index i.zip3
predict pg_resid, residual
rdrobust pg_resid running, h(37) all
scalar bpg_cct_2nd_1=_b[Conventional]
scalar sepg_cct_2nd_1=_se[Conventional]
rdrobust pg_resid running, h(52) all
scalar bpg_ik_2nd_1=_b[Conventional]
scalar sepg_ik_2nd_1=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pg_resid running, bwselect(CCT) all
rdrobust pg_resid running, bwselect(IK) all
preserve
gen zip2=(zip3-mod(zip3,10))/10
collapse pg_resid, by(running zip2 brand_index)
rdrobust pg_resid running, h(31) all
scalar bpg_cct_2nd_2=_b[Conventional]
scalar sepg_cct_2nd_2=_se[Conventional]
rdrobust pg_resid running, h(53) all
scalar bpg_ik_2nd_2=_b[Conventional]
scalar sepg_ik_2nd_2=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pg_resid running, bwselect(CCT) all
rdrobust pg_resid running, bwselect(IK) all
restore
* Results for Table A.6, price of E20/E25 gasoline, second discontinuity
dis "Coef: " bpg_cct_2nd_1 ", SE: " sepg_cct_2nd_1 
dis "Coef: " bpg_cct_2nd_2 ", SE: " sepg_cct_2nd_2 
dis "Coef: " bpg_ik_2nd_1  ", SE: " sepg_ik_2nd_1 
dis "Coef: " bpg_ik_2nd_2  ", SE: " sepg_ik_2nd_2 

use $fuelprices, clear
sort brand
gen brand_index=_n if brand~=brand[_n-1]
replace brand_index=brand_index[_n-1] if brand==brand[_n-1]
sort date_surveyed
gen running=date_surveyed-date("01-Oct-2011", "DMY")
keep if running>=-90&running<=89
drop if zip3==. 
quietly reg pg c.running##c.running##c.running##c.running##c.running##c.running##c.running##c.running ///
	   i.brand_index i.zip3
predict pg_resid, residual
rdrobust pg_resid running, h(32) all
scalar bpg_cct_3rd_1=-_b[Conventional]
scalar sepg_cct_3rd_1=_se[Conventional]
rdrobust pg_resid running, h(52) all
scalar bpg_ik_3rd_1=-_b[Conventional]
scalar sepg_ik_3rd_1=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pg_resid running, bwselect(CCT) all
rdrobust pg_resid running, bwselect(IK) all
preserve
gen zip2=(zip3-mod(zip3,10))/10
collapse pg_resid, by(running zip2 brand_index)
rdrobust pg_resid running, h(24) all
scalar bpg_cct_3rd_2=-_b[Conventional]
scalar sepg_cct_3rd_2=_se[Conventional]
rdrobust pg_resid running, h(69) all
scalar bpg_ik_3rd_2=-_b[Conventional]
scalar sepg_ik_3rd_2=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pg_resid running, bwselect(CCT) all
rdrobust pg_resid running, bwselect(IK) all
restore
* Results for Table A.6, price of E20/E25 gasoline, third discontinuity
dis "Coef: " bpg_cct_3rd_1 ", SE: " sepg_cct_3rd_1 
dis "Coef: " bpg_cct_3rd_2 ", SE: " sepg_cct_3rd_2 
dis "Coef: " bpg_ik_3rd_1  ", SE: " sepg_ik_3rd_1 
dis "Coef: " bpg_ik_3rd_2  ", SE: " sepg_ik_3rd_2 

use $fuelprices, clear
sort brand
gen brand_index=_n if brand~=brand[_n-1]
replace brand_index=brand_index[_n-1] if brand==brand[_n-1]
sort date_surveyed
gen running=date_surveyed-date("01-Feb-2010", "DMY")
keep if running>=-90&running<=89
drop if zip3==. 
quietly reg pe c.running##c.running##c.running##c.running##c.running##c.running##c.running##c.running ///
	   i.brand_index i.zip3
predict pe_resid, residual
rdrobust pe_resid running, h(10) all
scalar bpe_cct_1st_1=-_b[Conventional]
scalar sepe_cct_1st_1=_se[Conventional]
rdrobust pe_resid running, h(24) all
scalar bpe_ik_1st_1=-_b[Conventional]
scalar sepe_ik_1st_1=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pe_resid running, bwselect(CCT) all
rdrobust pe_resid running, bwselect(IK) all
preserve
gen zip2=(zip3-mod(zip3,10))/10
collapse pe_resid, by(running zip2 brand_index)
rdrobust pe_resid running, h(11) all 
scalar bpe_cct_1st_2=-_b[Conventional]
scalar sepe_cct_1st_2=_se[Conventional]
rdrobust pe_resid running, h(23) all
scalar bpe_ik_1st_2=-_b[Conventional]
scalar sepe_ik_1st_2=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pe_resid running, bwselect(CCT) all
rdrobust pe_resid running, bwselect(IK) all
restore
* Results for Table A.6, price of E100 ethanol, first discontinuity
dis "Coef: " bpe_cct_1st_1 ", SE: " sepe_cct_1st_1 
dis "Coef: " bpe_cct_1st_2 ", SE: " sepe_cct_1st_2 
dis "Coef: " bpe_ik_1st_1  ", SE: " sepe_ik_1st_1 
dis "Coef: " bpe_ik_1st_2  ", SE: " sepe_ik_1st_2 

use $fuelprices, clear
sort brand
gen brand_index=_n if brand~=brand[_n-1]
replace brand_index=brand_index[_n-1] if brand==brand[_n-1]
sort date_surveyed
gen running=date_surveyed-date("01-May-2010", "DMY")
keep if running>=-90&running<=89
drop if zip3==. 
quietly reg pe c.running##c.running##c.running##c.running##c.running##c.running##c.running##c.running ///
	   i.brand_index i.zip3
predict pe_resid, residual
rdrobust pe_resid running, h(10) all
scalar bpe_cct_2nd_1=_b[Conventional]
scalar sepe_cct_2nd_1=_se[Conventional]
rdrobust pe_resid running, h(21) all
scalar bpe_ik_2nd_1=_b[Conventional]
scalar sepe_ik_2nd_1=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pe_resid running, bwselect(CCT) all
rdrobust pe_resid running, bwselect(IK) all
preserve
gen zip2=(zip3-mod(zip3,10))/10
collapse  pe_resid, by(running zip2 brand_index)
rdrobust pe_resid running, h(11) all
scalar bpe_cct_2nd_2=_b[Conventional]
scalar sepe_cct_2nd_2=_se[Conventional]
rdrobust pe_resid running, h(21) all
scalar bpe_ik_2nd_2=_b[Conventional]
scalar sepe_ik_2nd_2=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pe_resid running, bwselect(CCT) all
rdrobust pe_resid running, bwselect(IK) all
restore
* Results for Table A.6, price of E100 ethanol, second discontinuity
dis "Coef: " bpe_cct_2nd_1 ", SE: " sepe_cct_2nd_1 
dis "Coef: " bpe_cct_2nd_2 ", SE: " sepe_cct_2nd_2 
dis "Coef: " bpe_ik_2nd_1  ", SE: " sepe_ik_2nd_1 
dis "Coef: " bpe_ik_2nd_2  ", SE: " sepe_ik_2nd_2 

use $fuelprices, clear
sort brand
gen brand_index=_n if brand~=brand[_n-1]
replace brand_index=brand_index[_n-1] if brand==brand[_n-1]
sort date_surveyed
gen running=date_surveyed-date("01-Oct-2011", "DMY")
keep if running>=-90&running<=89
drop if zip3==. 
quietly reg pe c.running##c.running##c.running##c.running##c.running##c.running##c.running##c.running ///
	   i.brand_index i.zip3
predict pe_resid, residual
rdrobust pe_resid running, h(21) all
scalar bpe_cct_3rd_1=-_b[Conventional]
scalar sepe_cct_3rd_1=_se[Conventional]
rdrobust pe_resid running, h(22) all
scalar bpe_ik_3rd_1=-_b[Conventional]
scalar sepe_ik_3rd_1=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pe_resid running, bwselect(CCT) all
rdrobust pe_resid running, bwselect(IK) all
preserve
gen zip2=(zip3-mod(zip3,10))/10
collapse pe_resid, by(running zip2 brand_index)
rdrobust pe_resid running, h(19) all
scalar bpe_cct_3rd_2=-_b[Conventional]
scalar sepe_cct_3rd_2=_se[Conventional]
rdrobust pe_resid running, h(27) all
scalar bpe_ik_3rd_2=-_b[Conventional]
scalar sepe_ik_3rd_2=_se[Conventional]
* CCT and IK bandwidth criteria
rdrobust pe_resid running, bwselect(CCT) all
rdrobust pe_resid running, bwselect(IK) all
restore
* Results for Table A.6, price of E100 ethanol, third discontinuity
dis "Coef: " bpe_cct_3rd_1 ", SE: " sepe_cct_3rd_1 
dis "Coef: " bpe_cct_3rd_2 ", SE: " sepe_cct_3rd_2 
dis "Coef: " bpe_ik_3rd_1  ", SE: " sepe_ik_3rd_1 
dis "Coef: " bpe_ik_3rd_2  ", SE: " sepe_ik_3rd_2 

end


************************************************************************************
************************************************************************************
************************************************************************************
capture program drop figures
program define figures
************************************************************************************
************************************************************************************
************************************************************************************

* This part reproduces Figure 4
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
gen time_1=-1-time
* Figure 4a
twoway lfitci res_o3_1st time_1 if time_1>=-30&time_1<0 || lfitci res_o3_1st time_1 if time_1>=0&time_1<=29 ///
	|| scatter res_o3_1st time_1 if time_1>=-30&time_1<=29 , msize(small) mstyle(p1) /// 
	legend(off) ///
	xlabel(-30 0 29, labsize(small) angle(0)) xline(0) ///
	ylabel(-2 -1 0 1, labsize(small) angle(0)) ///
	ytitle("log ozone residuals", size(medsmall)) ///
	xtitle("Days from E20 to E25 shift", size(medsmall)) scheme(s1manual)
* Figure 4b
twoway lfitci res_o3_2nd time if time>=-30&time<0 || lfitci res_o3_2nd time if time>=0&time<=29,  ///
	|| scatter res_o3_2nd time if time>=-30&time<=29 , msize(small) mstyle(p1) /// 
	legend(off) ///
	xlabel(-30 0 29, labsize(small) angle(0)) xline(0) ///
	ylabel(-2 -1 0 1, labsize(small) angle(0)) ///
	ytitle("log ozone residuals", size(medsmall)) ///
	xtitle("Days from E20 to E25 shift", size(medsmall)) scheme(s1manual)
* Figure 4c	
twoway lfitci res_o3_3rd time_1 if time_1>=-30&time_1<0 || lfitci res_o3_3rd time_1 if time_1>=0&time_1<=29,  ///
	|| scatter res_o3_3rd time_1 if time>=-30&time<=29 , msize(small) mstyle(p1) /// 
	legend(off) ///
	xlabel(-30 0 29, labsize(small) angle(0)) xline(0) ///
	ylabel(-2 -1 0 1, labsize(small) angle(0)) ///
	ytitle("log ozone residuals", size(medsmall)) ///
	xtitle("Days from E20 to E25 shift", size(medsmall)) scheme(s1manual) 
* Figure 4d
twoway lfitci res_o3_4th time if time>=-30&time<0 || lfitci res_o3_4th time if time>=0&time<=29, ///
	|| scatter res_o3_4th time if time>=-30&time<=29 , msize(small) mstyle(p1) /// 
	legend(off) ///
	xlabel(-30 0 29, labsize(small) angle(0)) xline(0) ///
	ylabel(-2 -1 0 1, labsize(small) angle(0)) ///
	ytitle("log ozone residuals", size(medsmall)) ///
	xtitle("Days from E20 to E25 shift", size(medsmall)) scheme(s1manual)  
	
*************************************************************************************
*************************************************************************************	

* This part reproduces Figure 5
* File bandwidth_sensitivity.dta must first be generated from the code below, 
* following "code for sensitivity analysis to variation in the optimal bandwidth"
/*
use bandwidth_sensitivity, clear
twoway line o3 o3_low o3_high bd if dis==1 & bd>=20 &bd<=41, ///
	   legend(off) lpattern(solid solid solid) lcolor(b b b) lwidth(medthick medthick medthick) ///
	   xlabel(20 41, labsize(small)) xmlabel(30 "30/CCT" 31 "31/IK", angle(270) labsize(small)) ///
	   ylabel(-0.15(0.05)0.25, labsize(small) angle(horizontal)) ///
	   xline(30 31, lpattern(shortdash) lcolor(b)) ///
	   yline(0, lpatter(shortdash) lcolor(b)) ///
	   ytitle("95% CI of treatment effect", size(medsmall)) ///
	   xtitle("Bandwidth in days", size(medsmall)) ///
	   scheme(s1manual)
twoway line o3 o3_low o3_high bd if dis==2 & bd>=15 &bd<=45, ///
	   legend(off) lpattern(solid solid solid) lcolor(b b b) lwidth(medthick medthick medthick) ///
	   xlabel(15 45, labsize(small)) xmlabel(25 "25/CCT" 35 "35/IK", angle(270) labsize(small)) ///
	   ylabel(-0.15(0.05)0.25, labsize(small) angle(horizontal)) ///
	   xline(25 35, lpattern(shortdash) lcolor(b)) ///
	   yline(0, lpatter(shortdash) lcolor(b)) ///
	   ytitle("95% CI of treatment effect", size(medsmall)) ///
	   xtitle("Bandwidth in days", size(medsmall)) ///
	   scheme(s1manual)
twoway line o3 o3_low o3_high bd if dis==3 & bd>=12 &bd<=34, ///
	   legend(off) lpattern(solid solid solid) lcolor(b b b) lwidth(medthick medthick medthick) ///
	   xlabel(12 34, labsize(small)) xmlabel(22 "22/CCT" 24 "24/IK", angle(270) labsize(small)) ///
	   ylabel(-0.15(0.05)0.25, labsize(small) angle(horizontal)) ///
	   xline(22 24, lpattern(shortdash) lcolor(b)) ///
	   yline(0, lpatter(shortdash) lcolor(b)) ///
	   ytitle("95% CI of treatment effect", size(medsmall)) ///
	   xtitle("Bandwidth in days", size(medsmall)) ///
	   scheme(s1manual)
twoway line o3 o3_low o3_high bd if dis==4 & bd>=13 &bd<=38, ///
	   legend(off) lpattern(solid solid solid) lcolor(b b b) lwidth(medthick medthick medthick) ///
	   xlabel(13 38, labsize(small)) xmlabel(28 "28/CCT" 23 "23/IK", angle(270) labsize(small)) ///
	   ylabel(-0.1(0.05)0.15, labsize(small) angle(horizontal)) ///
	   xline(23 28, lpattern(shortdash) lcolor(b)) ///
	   yline(0, lpatter(shortdash) lcolor(b)) ///
	   ytitle("95% CI of treatment effect", size(medsmall)) ///
	   xtitle("Bandwidth in days", size(medsmall)) ///
	   scheme(s1manual)

* code for sensitivity analysis to variation in the optimal bandwidth
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
collapse res_o3* time, by(date siteid dis)
forvalues var1=20(1)41{
	rdrobust res_o3_1st time, h(`var1')
}
forvalues var1=15(1)45{
	rdrobust res_o3_2nd time, h(`var1')
}
forvalues var1=12(1)34{
	rdrobust res_o3_3rd time, h(`var1')
}
forvalues var1=13(1)38{
	rdrobust res_o3_4th time, h(`var1')
}
*/

*************************************************************************************
*************************************************************************************	   

* This part reproduces Figure 6
* File demean_alternatives.dta must first be generated from the code below, 
* following "code for alternative procedures to demean the data in the first step"
/*
use demean_alternatives, clear
twoway line o3 o3_1 o3_2 bd if bd>=20&bd<=40&dis==1, legend(off) ///
	   lpattern(solid dash shortdash_dot) lcolor(b b b) ///
	   xlabel(20 30 40, labsize(small)) ///
	   ylabel(0 0.05 0.1 0.15, labsize(small) angle(0)) ///
	   ytitle("Treatment effect", size(medsmall)) ///
	   xtitle("Bandwidth in days", size(medsmall)) ///
	   scheme(s1manual)  
twoway line o3 o3_1 o3_2 bd if bd>=20&bd<=40&dis==2, legend(off) ///
	   lpattern(solid dash shortdash_dot) lcolor(b b b) ///
	   xlabel(20 30 40, labsize(small)) ///
	   ylabel(0 0.05 0.1 0.15, labsize(small) angle(0)) ///
	   ytitle("Treatment effect", size(medsmall)) ///
	   xtitle("Bandwidth in days", size(medsmall)) ///
	   scheme(s1manual)  	   
twoway line o3 o3_1 o3_2 bd if bd>=20&bd<=40&dis==3, legend(off) ///
	   lpattern(solid dash shortdash_dot) lcolor(b b b) ///
	   xlabel(20 30 40, labsize(small)) ///
	   ylabel(0 0.05 0.1 0.15, labsize(small) angle(0)) ///
	   ytitle("Treatment effect", size(medsmall)) ///
	   xtitle("Bandwidth in days", size(medsmall)) ///
	   scheme(s1manual)	   
twoway line o3 o3_1 o3_2 bd if bd>=20&bd<=40&dis==4, legend(off) ///
	   lpattern(solid dash shortdash_dot) lcolor(b b b) ///
	   xlabel(20 30 40, labsize(small)) ///
	   ylabel(0 0.05 0.1 0.15, labsize(small) angle(0)) ///
	   ytitle("Treatment effect", size(medsmall)) ///
	   xtitle("Bandwidth in days", size(medsmall)) ///
	   scheme(s1manual)

* code for alternative procedures to demean the data in the first step
* estimates directly on demeaned log concentrations at the site-hour-date level 
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
forvalues var1=20(1)40{
	rdrobust res_o3_1st time, h(`var1')
}

forvalues var1=20(1)40{
	rdrobust res_o3_2nd time, h(`var1')
}

forvalues var1=20(1)40{
	rdrobust res_o3_3rd time, h(`var1')
}

forvalues var1=20(1)40{
	rdrobust res_o3_4th time, h(`var1')
}	

* log concentrations are averages across hours within site-date
use $pollutants, clear
prepareVariables
quietly{
keep if hour>=12 & hour<=16
collapse (mean) O3 time tpl0 wsl0 rdl0 hml0 dvti0to199 dvti200to499 ppl0to3 ///
				kmcity_am pe_pg_ratio dayofweek weekday_reg ///
				wdqdt1l0 wdqdt2l0 wdqdt3l0 wdqdt4l0 beltway, ///
				by (date siteid dis)
gen logO3=log(O3)
foreach var1 in tpl0 wsl0 rdl0 hml0 {
	gen log`var1'=log(`var1')
}
gen logkmcity_am=log(kmcity_am+1)
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
			      c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&dis==1
predict res_temp if (siteid==`var1')&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&dis==2
predict res_temp if (siteid==`var1')&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&dis==3				  
predict res_temp if (siteid==`var1')&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time ///
				  c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio /// 
				  i.dayofweek ///
				  if (siteid==`var1')&dis==4	
predict res_temp if (siteid==`var1')&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&dis==4
drop res_temp
}
}
forvalues var1=20(1)40{
	rdrobust res_o3_1st time, h(`var1')
}

forvalues var1=20(1)40{
	rdrobust res_o3_2nd time, h(`var1')
}

forvalues var1=20(1)40{
	rdrobust res_o3_3rd time, h(`var1')
}

forvalues var1=20(1)40{
	rdrobust res_o3_4th time, h(`var1')
}
*/

*************************************************************************************
*************************************************************************************

* This part reproduces Figure 7
* File site_specific_regression.dta must first be generated from the code below, 
* following "code to implement second-step regression separately by site"
/*
use site_specific_regression, clear
hist Estimates if method=="CCT", xline(0, lpattern(dash)) start(-0.2) width(0.05) normal fcolor(none) lcolor(black) ///
	   addplot(hist Estimates if method=="CCT", start(-0.2) width(0.05) percent yaxis(2) ytitle("Percent", size(medsmall) axis(2)) fcolor(none) lcolor(black)) ///
	   scheme(s1manual) legend(off) ///
	   xtitle("E20 to E25 treatment effect, by ozone monitor-discontinuity pair", size(medsmall)) ///
	   ytitle("Density", size(medsmall)) ///
	   ylabel(0 2 4 6, axis(1) angle(0) labsize(small)) ylabel(0 10 20 30, axis(2) angle(0) labsize(small)) ///
	   xlabel(-0.2 -0.1 0 0.1 0.2 0.3, labsize(small))

hist Estimates if method=="IK", xline(0, lpattern(dash)) start(-0.2) width(0.05) normal fcolor(none) lcolor(black) ///
	   addplot(hist Estimates if method=="IK", start(-0.2) width(0.05) percent yaxis(2)  ytitle("Percent", size(medsmall) axis(2)) fcolor(none) lcolor(black)) ///
	   scheme(s1manual) legend(off) ///
	   xtitle("E20 to E25 treatment effect, by ozone monitor-discontinuity pair", size(medsmall)) ///
	   ytitle("Density", size(medsmall)) ///
	   ylabel(0 2 4 6 8, axis(1) angle(0) labsize(small)) ylabel(0 10 20 30 40, axis(2) angle(0) labsize(small)) ///
	   xlabel(-0.2 -0.1 0 0.1 0.2 0.3, labsize(small))

* code to implement second-step regression separately by site
use $pollutants, clear
prepareVariables
quietly{
gen res_o3_1st=.
gen res_o3_2nd=.
gen res_o3_3rd=.
gen res_o3_4th=.

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
			     (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==1
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1, residual
replace res_o3_1st=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==1
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg c.beltway ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==2
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2, residual
replace res_o3_2nd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==2
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==3				  
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3, residual
replace res_o3_3rd=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==3
drop res_temp
}
foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
reg logO3 (c.time##c.time##c.time##c.time##c.time##c.time##c.time##c.time)#i.hour ///
				  (c.logtpl0##c.logtpl0 c.logrdl0##c.logrdl0 c.loghml0##c.loghml0 c.logwsl0##c.logwsl0 ///
				  c.wdqdt1l0 c.wdqdt2l0 c.wdqdt3l0 c.wdqdt4l0 c.ppl0to3 ///
				  c.dvti0to199 c.dvti200to499 ///
				  c.logkmcity_am#c.weekday_reg ///
				  c.pe_pg_ratio##c.pe_pg_ratio##c.pe_pg_ratio)#i.hour /// 
				  i.dayofweek##i.hour ///
				  if (siteid==`var1')&hour>=12&hour<=16&dis==4	
predict res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4, residual
replace res_o3_4th=res_temp if (siteid==`var1')&hour>=12&hour<=16&dis==4
drop res_temp
}
}
foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
		rdrobust res_o3_1st time if siteid==`var1', bwselect(CCT)
		rdrobust res_o3_1st time if siteid==`var1', bwselect(IK)
}

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
		rdrobust res_o3_2nd time if siteid==`var1', bwselect(CCT)
		rdrobust res_o3_2nd time if siteid==`var1', bwselect(IK)
}

foreach var1 in 1 2 3 5 6 7 15 18 22 27 29 31 {
		rdrobust res_o3_3rd time if siteid==`var1', bwselect(CCT)
		rdrobust res_o3_3rd time if siteid==`var1', bwselect(IK)
}

foreach var1 in 1 2 3 5 6 7 15 16 18 22 27 28 29 31 33 35 37 {
		rdrobust res_o3_4th time if siteid==`var1', bwselect(CCT)
		rdrobust res_o3_4th time if siteid==`var1', bwselect(IK)		
}
*/

end
